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priate to  the  processing  of  seismic  data.  One  is  its  ability  to  apply  a given 
algorithm  to  sixty-four  different  data  streams  simultaneously,  thus  providing 
an  order  of  magnitude  increase  in  processing  speed  over  conventional  machines. 
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ABSTRACT 


Two  features  of  the  ILLIAC  [V  system  at  NASA/Ames  are  particularly  appro- 
priate to  the  processing  of  seismic  data.  One  Is  Its  ability  to  apply  a given 
algorithm  to  sixty-four  different  data  streams  simultaneously,  thus  providing 
an  order  of  magnitude  increase  in  processing  speed  over  conventional  machines. 
The  second  is  its  large  data  storage.  The  seismological  algorithms  for 
convolution  filtering,  beamforming,  matched  filtering,  PHILTRE,  maximum  likeli- 
hood, and  FKCOMB  are  each  able  to  take  advantage  of  these  features  in  the 
processing  of  seismic  data.  The  data  can  be  arranged  so  that  each  processing 
element  contains  successive  time  windows  on  a given  trace,  as  in  bandpass 
filtering,  or  successive  beam  seta,  as  in  beamforming. 

Some  preliminary  data  editing  is  required  for  each  of  these  algorithms 
to  arrange  the  data  appropriately  in  processing  element  memory  to  utilize  the 
ILLIAC  IV  computer  efficiently.  Data  formatting  schemes  were  designed  for 
one  algorithm  (FKCOMB)  which  was  coded  and  implemented  on  the  ILLIAC;  these 

schemes  can  be  appropriately  modified  for  use  with  other  seismological  algo- 
rithms. 

Experience  with  data  transfer,  program  entry  and  editing,  compilation, 
and  program  execution  show  that  while  the  ILLIAC  system  is  still  under 
development,  adequate  facilities  do  exist  for  development  of  seismological 
algorithms. 
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INTRODUCTION 


The  purpose  ol  this  study  was  to  determin-i  the  suitability  of  the  ILLIAC 
computer  for  processing  seismic  data.  We  have  done  this  by  looking  at  the 
computing  requirements  of  each  of  several  algorithms;  and  then,  by  comparing 
these  requirements  with  the  characteristics  of  the  ILLIAC,  we  have  investigated 
the  feasibility  of  programming  each  of  the  algorithms  on  the  ILLIAC.  Finally, 
the  procedure  FKCOMB  was  actually  coded  for  the  ILLIAC  and  the  program  has 
been  tested  and  run.  FKCOMB  is  a long-period  seismic  signal  analysis  p-oce- 
dure,  which  is  important  in  calculating  discriminants  between  earthquakes  and 
nuclear  explosions;  it  may  become  an  integral  part  of  data  processing  on  the 
Seismic  Network.  FKCOMB  was  chosen  for  this  experiment  because  the  large 
amount  of  processor  time  required  prohibits  its  use  in-house.  Also,  known 
results  are  available  with  which  to  compare  the  ILLIAC  version. 

The  ILLIAC  computer  consists  of  a control  unit,  64  arithmetic  units  or 

9 

processing  elements  (PE),  128K  64-bit  words  of  core,  and  10  bits  of  disk 
storage.  The  64  PE's  execute  instructions  in  lock  step;  i.e.,  they  all 
execute  the  same  instruction  simultaneously.  It  is  in  this  respect  that  the 
ILLIAC  departs  from  conventional  computer  architecture. 

The  control  unit  decodes  Instructions  and  executes  instructions  for  pro- 
gram control.  It  has  24-bit  Integer  arithmetic  hardware  to  calculate  indices 
and  addresses.  There  are  four  general  purpose  accumulators  and  a 64-word 
fast  access  memory  in  the  CU,  which  also  has  access  to  all  128K  of  core  and 
initiates  transfers  between  core  and  disk. 

The  primary  computational  resource  of  the  ILLIAC  is  the  array  of  64 
processing  elements.  Each  has  complete  arithmetic  capabilities  and  can 
perform  2 x 10^  multlpllcetlons  per  second.  The  capacity  of  all  64  PE's  is 
about  10®  nultiplicatlons  per  second.  Each  PE  has  access  only  to  2K  of  core,  . 
and  has  only  limited  capability  to  communicate  with  other  PE's.  Control 
within  a PE  consists  of  the  ability  programatically  to  disable  selected  PE's. 
When  disabled,  a PE's  memory  is  protected  and  cannot  be  altered  by  the  PE, 
though  all  other  facets  of  instruction  are  performed. 
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Ilu*  II.I.IAI  ilisK  i.s  tliL*  prin.irv  .'itor.u'i'  licvici*.  It  consists  of  11  Wi'.ul- 

pc  r-t  riU.l;  iiisks  ainl  two  disk  cont  rol  Uts.  I opei  iior  tho  disks  holil  approxi- 

ir..itolv  10  bits.  Iransters  botwcun  toie  ami  di^>k  .ire  initiated  by  the  IldlAC 

v'ontrol  unit  and  cur  ui  bioi  ks  or  pap,»“s  t)l  i(Jd4  worils.  I'he  transler  rate 

>) 

is  about  111  bits  per  second.  Access  time  to  .in%  record  on  disk  is  AU  milli- 
seconds or  less.  Itie  disk  can  be  loaded  ! rom  the  lenex  !ile  svstem  prior  t',' 

pro);ram  execution  and  unloaded  alter  piopra-.  termination.  llie  layout  cjl  data 
on  the  IM.IAI  di-.k  is  under  user  corttrol  .iiui  should  be*  arranpc*d  so  as  to 
Tiininii'e  access  limes  duriny  proyrar  e xecut  i i)ii . 

idle  h4  processiiu’  element:,  provide  tiie  .■liiilii',  to  pertarm  vector  arith- 
••letic  operations  on  iSA  data  elerents  ; imu  1 1 .ineo*:;,  1 . Provran  loylc  generally 
reciuiros  that  selc'cted  1’1's  he  cHsableci  to  .c . lid  re*lundant  calculations  ji  all 
h’>  iirocessiin*  units  are  not  reciuired.  In  comral,  program  execution  time  is 
decreased  if  disabling  ot  i'K's  is  a*.'oick-*d. 

ILI.IAl'  is  able  to  periorr  aiiprox ic.atc  1 v I'lO  million  operations  Ij.e.,  a 
multiplication,  addition,  etc.i  per  second.  Any  proceclure  which  rec^uiies 
lewer  than  ino  tillion  operations  would  tiave  a runninp  time  ol  under  10  ninutes. 
idle  setup  time  lor  an  lld.IAl  job  is  lar>’e  enou.  b to  make  such  a run  impractical, 
dhus,  al.'or  i th'-,s  requirini*  verv  tew  computations  or  the  use  of  a small  data 
base  with  anv  .ilyorithm  are  unsuitahlc*  lor  lld.lAb  processing. 


AlM’KDACll  !0  I)i;Sl(.NlN(.  I’AKALi.KL  SKI  SMOLOC ICAI.  AI.t.OKITHMS 


Ovfrvlew  ol  tl.i.lAC 

II.I.IAt;  i.s  .1  parallel  processor.  It  consists  of  a control  unit  (CL'), 

0-4  proi  essin>;  elements  (PK)  , 1 31,072  words  of  core  memory,  and  15,974,^00 
words  of  dls'it  racmorv.  The  control  unit  has  access  to  all  of  core  memory. 

Us  basic  cycle  time  Is  bO  nanoseconds.  However,  greater  processing  power 
Is  achieved  through  tie  simultaneous  execution  of  an  instruction  In  each  of 
the  n4  processing  elements. 

The  control  unit  fetches  and  decodes  all  instructions.  After  decoding, 
some  Instructions  are  broadcast  for  execution  In  the  processing  elements  while 
otliers  are  exet  ited  in  the  control  unit.  The  arithmetic  capability  of  the 
control  unit  is  limited  to  24  bit  two's  complement  addition  and  subtraction, 
masking,  and  com;iarison  for  use  in  branching.  The  control  unit  has  no  floating 
point  .apability.  One  operand  at  .a  time  is  processed  by  the  control  unit. 

The  control  uni:  also  initiates  data  transfers  between  core  and  ILl.lAC  disk. 

The  processing  power  of  ll.l.iAC  resides  in  b4  identical  processing 
elements.  r.ach  I’h  ext^cutes  Instructions  broadcast  from  the  CL.  Though  each 
I’K  has  its  own  index  registers  and  memory  to  operate  upon,  all  64  PE's  always 
execute  identical  instructions  in  lockstep.  Each  PE  has  direct  access  to 
2048  words  ol  core  memorv,  shown  as  a column  in  Figure  1. 

There  are  three  data  paths  available  for  communication  among  PE's  and 
between  the  I’E's  and  the  control  unit  (CL).  First,  the  CL  can  access  all  of 
core,  so  It  can  load  a word  from  one  processing  element  memory  (PEM)  and 
either  use  It  or  store  it  in  another  PEM.  This  method  of  communication  is 
hcJth  simple  and  flexible,  allowing  for  any  data  movement  desired,  but,  since 
only  one  word  at  a time  is  transferred,  it  is  relatively  slow  compared  to 
the  two  other  methods  available. 

Second,  the  (T  can  communicate  with  all  PE's  by  broadcasting  the  same 
word  to  all  I’h's  simultaneously.  This  method  is  faster  than  the  first  since 
sixtv-four  words  are  transmitted  at  once,  but  provides  only  a limited  form  of 
common  1 cat  ion. 
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Hurd,  the  PE's  can  con»u„lcate  with  each  other  via  the  ROUTE  Instruction 

i';oL:::n::chtr°":r,r:n:  ^ 

irr:  a-'r  — ::::::::: 

huniber  of  the  destination  PE  is  Moll  (PEK+R)  If 
Ehe  PE's  are  thought  of  as  arranged  In  a circle  with  PF  61  H “ 

--  - consists  Of  loading  the  data.  llZg^L  :::::  Id" 

«or  „g  the  data.  This  data  transfer  Is  very  fast  since  64  words  are  trans- 
erred  simultaneously.  It  Is  general  m that  all  64  words  can  he  different 
cue  pattern  set  hy  the  fact  that  the  routing  distance  Is  the  same  for  11 

cor  hTdlT'"' 

core  to  64  different  locations  simultaneously. 

The  primary  memory  used  by  ILLIAC  is  a disk  memory  with  capacity  approk- 

ima  e y 00  times  that  of  core  memory.  One  page  (1024  wo-d  i 
, UU24  words)  of  memory  is 

the  minimum  amount  of  data  that-  rrnr.  - n 

>>c  transferred  between  core  and  disk. 

Although  the  bandwidth  between  core  and  disk  is  5yin^  Ku 

UJ.SK.  IS  .5X1U  bits  per  second  thp 

average  access  time  to  a particular  spot  on  disk  is  20  mill • 

relativplv  l >-  on  disk  is  20  milliseconds.  This 

parallel  pr  “ “ "‘■”‘’*"0"''  clock  time  in  64 

Ota::::::"'  -»^cr 

transfers  must  be  kept  to  a minimum  to  avoid  waiting  for  disk  accesses. 

Since  the  most  important  feature  of  Ihiun  Is  its  computing  power 
( pproximately  100  million  multiplications  per  second),  one  of  the  prime 
ectlves  in  the  design  of  any  llUAC  program  Is  to  minlmUe  execu  Ion  time 
e best  possible  result  is  a running  time  1/64  that  possible  with  only  Ine 
b,  but  due  to  the  architecture  of  the  machine  the  degree  to  which  this  is 
leved  Is  dependent  upon  the  design  of  the  algorithm.  He  discuss  two 
general  approaches  to  designing  algorithms.  First,  suppose  that  it  Is 
necessary  to  code  the  trigonometric  SIN  function  for  lUIAC.  If  the  particular 
u age  makes  It  possible  to  always  compute  64  functions  simultaneously,  one 
simple  has  the  same  SIN  routine  running  In  all  PE's  on  different  data  and 
= speedup  by  a factor  of  64  Is  very  nearly  achieved.  (Some  111:  TLTu 
ere  is  conditional  branching  in  the  original  SIN  routine  which  Is  changed 
enabling  and  disabling  of  PE's).  A second  approach  is  to  devise  a method 
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or  uuUxans  64  ,,r..C6s6„r«  to  c.„,.,uc.  v„l,.c.  N„  ,„u,„l 

7,  ""■  •“*  “ O'  •«"  » fa«6r  tl„.„  „ 

with  „„6  nr.r  a,.,. reel,  b„c„  faata,-  a,.„Ur,  but 

curtatn  alsortttaa  „ay  ptacludo  calculation  of  ,»>rc  tl.an  one  value  of  SIN 
s.n,ulta„coualv  or  .ay  ro,„irc.  sl„„iflcant  overhead  ol.soohnro  In  order  to  do 

SO. 

one  .laconcuptlon  la  that  if  „1,  „f  the  Pp.a  ,,ro  Itopt  "hnay"  the  .nchine 
ta  run„.„8  at  «xd,„„  efficiency.  In  fact  this  .statement  is  not  true  and  one 
must  he  very  careful  relatl„«  the  cord  efficiency  to  the  use  of  imAC 
tor  example,  consider  the  prohle.  of  summing  groups  of  numbers.  If  it  is 
desired  to  .sum  64  pairs  of  numbers,  keeping  64  different  re.sults,  each  PF 
rorms  one  sum  and  the  work  Is  dene  64  times  faster  than  could  be  done  by  one 
processor.  if  however,  it  Is  desired  to  find  the  sum  of  on,  group  of  64 

numbers,  a more  complicated  method  must  be  u.sed.  In  order  to  simplify  the 

explanation  somewhat,  consider  an  oioKr  di.-  i • 

. consider  an  eight  PL  machine  and  the  summation  of  eight 

numbers,  one  in  each  PE.  Figure  2 depicts  a method  whereby  this  can  be  done 

in  three  routes  and  three  additions,  .since  the  routes  reguire  roughly  egul- 

V ent  CPU  time  as  the  register  loads  necessary  before  any  operation,  the 

time  taken  lor  .™  8 PE  machine  to  sum  eight  numbers  is  equal  to  the  time 

taken  for  three  additions,  if  this  algoriti™  is  extended  to  the  summation 

of  64  numbers  within  a 64  PE  machine  it  takes  6 additions  to  form  the  sum. 

■iven  that  one  PE  requires  63  .additions  to  sum  64  numbers,  the  64  PE  machine 

1»  63/6  or  10.5  times  faster.  .Note  that  although  none  of  the  PE's  are  ever 

disabled  and  all  .are  forming  tin.,  sum,  this  algorithm  does  not  achieve  the 

factor  of  64  speed  up.  however,  the  factor  of  ten  speed  up  that  is  achieved 

makes  this  algoriti™  usable  if  data  organisation  requires  its  use. 

Dc.sign  Considerations 

The  choice  of  which  design  approach  to  t.ake  for  a particular  problem  is 
ependent  upon  data  organisation.  There  is  often  one  approach  requiring  a 
very  specific  data  organization  which  Is  much  faster  than  any  other.  It  must 
be  decided  whether  the  overhead  and  execution  time  involved  In  data  trans- 

pos.tlon  Is  compensated  by  decreased  overhead  and  execution  time  elsewhere 
in  the  algoritlim. 
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STATE  OF  RF.(;iSTF.KS 
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M>*,iire  J,  ot  Ki’wsun  ion. 


Ttui  lirst  step  in  desi^;nin^!  a seismic  aigorithin  to  run  on  II.LIAC  is  to 
examine  similar  or  repeated  data  structures  and  determine  how  they  could  be 
organizeit  in  the  processor  memory  and  to  analyze  similar  or  repeated  opera- 
tions and  determine  tiow  they  could  be  divided  among  tlie  processors. 

I.ong  and  short  period  seismic  data  are  recorded  at  seismic  arrays  con- 
sisting of  a group  ol  sensors  samjjled  at  a constant  time  interval.  The  data 
so  recorded  consists  of  a series  of  data  scans.  bach  data  scan  is  a time 
sample  1 rom  each  sensor.  There  are  two  sti'iu  tores  repeatetl  througlujut  the 
data.  lirst,  there  are  several  channels,  each  identical  in  structure  to  the 
rest.  Second,  there  are  many  identically  structured  time  samples.  in  order 
to  utilize  either  of  these  structures,  time  must  be  spent  transposing  the 
data.  It  would  bo  convenient  it  it  were  possible  to  process  the  data  without 
transposing  in  any  way  - but  the  input  consists  of  data  records  which  are 
formatted  differently  for  each  array. 

Since  the  data  must  be  restructured,  it  is  reasonable  to  build  a new 
structure  which  makes  processing  as  fast  and  straight  forward  as  possible. 

The  choice  between  the  two  data  structures  is  dependent  upon  the  requirements 
ol  the  algorithm.  General  discussions  ol  the  several  seismic  algorithms  and 
their  data  requirements  is  contained  in  Section  3.  A detailed  discussion  of 
the  design  ol  FKCOMb  is  found  in  Section 
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SFISMK;  AITI.lCA'l  ions  on  i'li!';  II.I.IAC  IV  CiiMPOTKi'' 


OcMiera  L 

In  (,1r‘  following’  pai  .'4.'.  (li.si.iiss  liu.‘  snitahilitv  ol  Liu- 

II.I.IAC  comput-i'r  for  pia-icossinp  .sfi.anic  data  iisinj.;  .sovi' r.i  1 Icsinl  a l>.;or  i thins  . 

tXir  i nvost  i 4at  i on  has  rovoali  d Lliat  l hi*  II.I.IAC  coinpntcr  is  p.i'fimal  1> 
.siiitotl  for  procossiii);  of  soismir  dai.i  which  involves  frequently  ri'insa  1 ed  or 
slmnltaneous  identical  operations  using  dilierent  stt^;  ol  liata,  .and  can 
be  iirogramnied  in  .such  a way  that  the  proces.sinj;  is  [lertornied  s imu  1 1 aneous  Jy 
in  the  64  jn'oeess  ing  olerients  ol  this  coniputer.  Thus,  i!  desirable,  it  will 
lx‘  fea.sihlo  to  use  II.I.IAC  to  I'l'is'ess  routinely  all  long-iit  r iod  data  for  the 
lilannt'd  Seismic  Network.  In  addition  it  can  al.so  he  n.a-d  for  oil -line  pro- 
ce.ssing  of  selected  data. 

In  this  di.scussion  we  shall  concentrate  on  the  possibilities  ol  this 
computer  tor  the  dotestiun  and  discrimination  of  .seismic  events  usin>;  seismic 
array  dat.a.  Tiie  computer  can  also  he  used  in  other  seisinii  applications  too 
numerous  to  treat  iiere.  Seismic  arrays  record  the  earth  motion  in  two  sepa- 
rate iretpiency  band::,  slu't  t-|  t r i od  and  long-period,  wln'ih  lor  soiiie  purj'oses 
require  different  treatments  ix'cause  ol  the  ditlereiit  nature  ol  seismic  waves 
recorded  in  tiie  two  Itaiiiis.  .Some  ol  the  procacsses  discu.ssed  are  used  for  dat.i 
in  both  Ivinds  while  others  are  (ommonly  u.-.ed  only  for  It'nc  or  short  period 
d a t a . 

liie  most  common  signals  lor  I nvi  s t i c,a  t i or.  in  the  short  peritel  baiv.l  are 
ti.e  siiort-iier  iod  body  wave‘:,  particularly  the  shor  t-per  i lul  I’  lirsl  arriv.il. 

P waves  can  arrive  at  a seismic  st.it  ion  with  a wide  range  (U  apparent  '.'elo- 
cities  and  1 rom  all  I'ossihle  azimuths.  .Since  tie  bandwidth  ol  the  si);nal  is 

limited,  fretiuency  liltering  tends  to  enliance  tiie  .sign.il/noisi  r.itio.  The 

detect  ioi  threshold  in  tiie  short-per  U>(1  hand  is  low  relative  to  that  ol  the 
long-{ieriod  hand,  and  events  arc  usiiallv  detected  in  this  b.utd . 'I  he  .irrival 

azimuth  and  ajiiiartnt  veh'citv  0!  the  siior  t-j>er  i od  I’  waves  .it  an  array  yields 

a iirelirainary  eiiicenter  location  wiiich  can  lie  used  to  nari'ow  the  searcii  for 
waves  in  the  1 ong-jier  1 od  band.  In  the  1 ong-iu’r  iod  h.uid,  long-jie  r i od  body 
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waves  ami  surface  waves  are  the  si>;nals  of  interest.  When  used  in  conjunction 
with  short-period  data,  they  are  al]  proven  or  potential  discriminants  between 
explosions  and  earthquakes.  The  most  important  of  these  is  the  loriK-period 
KayleiRh  wave,  whicl  is  used  in  the  discriminant.  The  kayieiRh  wave 

has  several  characteristics  which  can  he  utilized  by  detection  alRorithms: 

1.  Waveform  (path-dependent); 

2.  Particle  motion; 

1.  Azimuth  and  apparent  velocity. 

Since  in  most  cases  detection  already  has  occurred  on  the  short-period  data, 
it  is  only  necessary  to  prove  or  disprove  the  presence  of  RayleiRh  or  other 
long-period  waves  arriving  from  roughly  the  direction  of  the  preliminary 
epicenters,  and  to  measure  the  wave  amplitude  if  [ircsent. 

The  following  seismic  processing  algorithms  will  he  discussed  in  this 

report : 

1.  Freipiency  (convolution/recursive)  filtering 

2.  ileamforming 

■J.  Matched  filtering 

4.  PMII.THE 

5.  Maximum  likelihood  f— k spectra 

6.  FKCOMli 

Ihe  last  four  of  these  have  only  limited  or  no  application  for  the 
short  period  hand.  Une  processor  (IKCOMH)  is  discussed  in  more  detail  since 
it  was  selected  to  be  demonstrated  on  the  1 1.1, 1 AC. 

Convolution  and  Recursive  Filtering 

Simple  filtering  is  the  convolution  of  a sei.smic  trace  with  some  arbitrary 
function  which  limits  the  bandwidth  of  the  output.  Recursive  filtering 
accomplishes  the  same  result,  hut  makes  use  of  a 'eedback  loop  to  reduce  the 
number  ot  arithmetic  operations  required. 

This  operation  can  he  represented  mathematically  in  the  form: 

m k 
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wliCM'e  all  indici's  aro  iruoi't-rs  wr  > t i • . 

^ ’ "i  ‘'iJiies  ot  tilt'  orij'iiril  cli);iLixod  trari' 

-d  y.  are  value.  ol  the  tiltered  output,  aud  a^.  h..  n.  and  h are  constants 

Me  cl.ou-e  01  which  is  dependent  on  the  tilter  lunction  to  he  perfonned. 

The  ll.UAi:  is  well  suited  to  perlorn  convolution  or  recursi  ;e  lilterinp 
simultaneously  on  all  processing  units.  Ihese  algorithms  can  he  used  lor 
liUerini;  all  elements  ol  an  artay  usi„,;  the  sa,„-  in.,,,  ,,, 

li^Tiited  si>.nal  in  wideband  noise,  or  utilised  i or  lilterinn  the  same  trace 
with  a set  01  lilters  to  perlonn  a last  1-ourier  analvsis  or  to  compute  spec- 
tral ratios  ior  d i sc  r ini  na  t io  , . he  parallel  ajporithm  tan  also  be  used  to 
simultaneously  deconvolve  si.t;.-iour  seismic  irate:,  remove  instrument 
response,  simulate  seismuprams  produced  hv  diilerent  instruments,  or  to 

reduce  the  seismopram  traces  s imul taneous 1 v to  accelerations,  velocities,  and 
displacements  ns  functions  of  time. 

h 1 is  a Hclic.natlc  i eprosentat  ion  ol  possible  i r ranpemen  ts  ol  data 
in  tl,«  II.IIAC  ™,-v  ,„r  Tc, aiu-n,,,;.  ,,  ,„f,erent 

channel  nr  data  i,  to  each  I'h,  with  the  »a;«  llltor  aidSicd  u,  all  I'tS. 

In  H,;urc  lb,  a Riven  channel  oi  darn  is  ln|nit  to  all  I'lt's,  with  a diilerent 
niter  applied  to  each  I'h.  Ilnnre  Ic  reprosoius  a coiiibinal  i on  of  the  previous 
esanples  ehich  the  I'h's  are  parti. ionej  i„t„  sev.ral  set.s,  ai,  the  I'f, 

in  a R.ven  set  receivinp  the  s..r,e  data  bat  operatli,,;  eith  diilerent 

1 liters. 

heaml  or:-, in^ 

.he.nniomln,;  Is  the  pr..,.„ss  t live- „ I , tin,;  sevetai  channels  ol  arrav 
data  and  sn.unln.i  thee  to  fore  a sinple  channel.  The  tloe  »1„  , ts  chosen  .are 
the  natural  delav,  in  tiee  ol  ...  rival  a hvpothetical  sl.inal  . rossi,,,;  tin 
array.  The  delays  are  lelined  win,  respei  , to  ,so«.  ...bitrarv  point  in  space, 
lor  pl.ine  waves  ol  constant  vein,  if, ■,  the  delavs  are 


, " r . 

1 i 


-diere  1 Is  the  Indcs  „i  the  i‘"  sensor,  r'  1,  the  location  ol  the  sensor  and 
the  slowness  of  tlio  sijpial  is: 

S = VV(V-Vj 
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wlioro  V is  the  velocity  ol  the  sifUK.I  across  tlie  arr.i'.  Ii  one  has  e„n,put.<I 
delavs  1 rom  the  true  S ol  a given  . Ignal,  tiiat  is,  1 rora  its  tine  speed  urui 
arrival  azinmtli,  and  lias  assimied  that  the  sig.nal  wavelorn  does  not  vary 
during  transit,  the  eiiect  ol  time  sliiiting  is  to  make  all  -he  h.uiuels  appear 
to  have  been  recorded  at  the  arbitral  v reference  point.  The  eiiect  ol  su:.mnng, 
therelore,  will  be  to  add  the  signal  to  itself  h- 1 times,  when.  X is  the  num- 
bei  ol  channels.  I he  signal  is  thus  reinl'orced.  It  the  noisi'  is 
random  and  uncorrolated  between  array  elements,  it  is  reduced  to  of  its 

oiigiiial  amplitude  by  the  -lummatirn.  Thus  heaml  orm  i n)>,  has  the  luntion  ol 
liK'tcMsinp  the  offoctivo  s i 1 -Lt)-no  i .‘-ro  r;itiii, 

One  can  estimate  the  speed  and  direction  of  propagation  ol  signaLs  bv 
landing  ihe^maximuni  of  the  time  average  of  the  squared  beam  values  rdenotod 
K")  on  the  a-|)laiie: 


■'-1  1 N 

1 lA 
i = 0 


i I 


d-1 


1 = 0 


where  B.  is  the  expression  in  bracket.,  the  bean,  oi  the  array;  x.  are  the 
i channel  data  samples;  t is  tlu'  sampling  interval,  j is  the  time  index; 

1 is  tlie  number  oi  time  points  over  whi.l.  ateraje  is  takm. 

Ihe  iirobabi  1 Itv  of  tlie  preseiue  of  the  signal  can  lie  estimated  by  the 
statistic  _ 

-■■A' 

■■a,-,.,' 

Where  the  denominator  is  the  time  average  of  the  .cin,  o-  He  square  (or  power) 
of  tlie  individual  input  traces  x,  after  the  beam  is  subtracted. 

Ibis  statistic  is  distributed  ripjiroxirnatel v as  a non-central  F variable 
with  degrees  oi  I reedom  determined  by  the  number  of  channels,  bandwidth,  an-1 

the  Lime  lengtli  of  averaging  (assuming  that  only  uncorrelated  noise  is  pre- 
sent) . 


Ihe  standard  i tables  can  he  usi'd  to  determine  the  s i gn  i 1 i caiice  of 
detection,  or  the  detection  can  he  automated  (Blandioid,  19/1). 
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iho  beams  can  also  be  displayed  tor  the  visual  detection  ol  the  waves 
ol  interest.  I-or  detection  oi  siirlace  waves  and  the  measurement  ol  M this 
i-s  .till  the  best  procedure.  Experienced  analysts  can  recoKui^.e  and  Lasure 
seismic  wave  arrivals  in  many  cases  where  automatic  machine  detection  schemes 
tail.  Routine  computation  ol  lone-period  beams  and  their  storage  in  the 
mass  store  event  files  wouid  be  a vaiuablc  routine  function  for  the  NHb  and 
would  require  a substantial  computational  power  easily  met  by  the  Il.I.lAC. 
Therefore  beamforminy  might  well  he  the  single  most  useful  algorithm  for 
implementation  on  the  M.LIAC  computer. 

Several  basic  computational  configurations  can  be  used  in  beam  processing. 
There  are  shown  in  Figure  A.  in  the  first  configuration  (Aa)  each  IF  contains 
one  sensor  trace  and  tl.e  beam  values  are  accumulated  by  forming  row  sums  on 
the  se.eral  ! I.  s.  This  configuration  is  suitable  to  process  several  arrays 
simultaneously,  and  computing  the  desired  beams  from  a single  data  set 
sequentially  may  use  long  time  windows  such  as  those  needed  for  the  recogni- 
tion of  dispersed  surface  wave  trains.  Another  advantage  of  this  configuration 
is  that  preprocessing  of  traces  such  as  filtering  processing  can  also  be 
performed  simultaneously  prior  to  beaming  without  the  need  to  remove  the  data 
from  memory.  The  output  of  such  a scheme  can  he  directly  used  in  network 
event  processing.  This  configuration,  uniquely  possible  on  the  ILLIAC  IV,  is 
the  most  efficient  if  the  maximum  number  of  FF's  can  be  utilized.  This  can 
be  achieved  if  the  total  number  of  sensors  in  the  arrays  are  close  to  sixty- 
lour,  or  aiternatively  the  remaining  PF’n  are  used  to  compute  different  beams 
on  the  same  arrays.  To  obtain  continuous  seismograms  of  long  duration  this 
seems  to  be  the  most  efficient  approach,  since  various  preprocessing  schemes, 
sucl  as  convolution  filtering,  coordinate  rotation  to  obtain  Rayleigh  hove, 

SH  and  SV  components  can  he  performed  on  them  without  the  need  for  excessive 
numbers  of  overlaps  in  the  successive  time  windows  which  would  be  required 
if,  as  we  discuss  below,  each  I'F  were  to  contain  all  the  channels  of  data 
required  lor  a particular  beam.  Incidentally,  PIIII.TRF  can  be  used  as  a 

postprocessor  for  64  array  beams  previously  obtained  (for  64  events)  which 
can  be  run  in  parallel. 
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rh.Mf  .III-  two  other  .ilternate  hut  .■ener.i ! 1 y les.s  eifeitive  ( oiiiputa- 
tional  con!  ij’.iirat  ions  which  are  imlicated  in  Fipure  U.  one  loads  all  sensor 
traces  iron  one  arrav  into  one  ih  ami  ea.  h I’h  contains  a diiterent  tine  window. 
The  desired  beans  for  a piven  tine  window  na’.  th.-n  he  coiiputed  s.a|uent  i al  ly 
{Figure  4h).  The  other  sthene  Ki/ure  /,  • ) loa  is  the  sane  tine  window,  all 
traces,  into  as  nany  i’l.'s  as  there  are  do  , i red  heans  and  the  beans  are  computed 
s i nu  1 1 .ineous  1 V . Ihe  d i sail  vati  t a.r,e  oi  ttu-  last  nesitioned  neiinids  is  that  since 
each  i’i  contains  ,ill  traces  the  cor  respond  i nr,  time  windows  must  he  shorter 
due  to  I’K  r;er.!orv  limitations.  itiis  pioiessin.  , iiuludinj;  he.nr.inr,,  will 
require  noi  e l onplicated  huiterinp  s^  iie.-.t  -,  !)etween  core  and  disk.  Therefore 
it  seems  tiiat  tlie  first  conputat  iuna  I scheme  lias  tin*  most  practiial  value, 
altliouqh  the  others  mav  be  use  i advant ayeous i y , tor  instani e,  for  enhancing 
short  hoilv  wave  arrivals.  The  maximum  utiliration  ot  the  lompiiter  requires 
the  consideration  oi  tlie  tvpo  oi  proiessinp  required,  numher  oi  traces  or 
arrays  and  the  length  of  tir.e  windows  to  be  ptoces-.ed. 

Matched  Filtering 

lliis  technique  utiiijies  tlie  wavelorm  ot  the  signal  to  be  detected 
(/liexander  and  kabenstine,  lhh7a,b),  Tlie  exjiected  waveform  oi  the  signal  is 
used  on  the  seismic  trace  as  a convolution  filter.  Ideally  the  expected 
wavelorm  is  identical  to  the  actual  one  and  in  the  resulting  tuiiput  trace 
t.ie  sii'.nal  is  irans’ormed  into  a pulse  wi.ich  is  shaped  lire  the  autocorrela- 
tion oi  the  sigiial  waveform.  In  pra  tlce  it  is  not  possiiile  to  predict  the 
actual  wavelorm  preciselv,  so  itie  matched  filtering  results  in  the  contraction 
Oi  the  actual  signal,  whiiti  for  surface  waves  can  fie  a long  wave  train,  into 
■a  mu(  h shorter  wavelorm.  by  compressing  tfie  same  amount  of  energy  into  a 
sfiorter  time  interval,  liie  signal/noise  is  increased.  ft  also  de-emphasi xes 
sic.nals  which  do  not  natch  the  wavelorm  used  :or  filtering.  Ihe  technique 
ii.is  been  su.cessiul  in  detecting  .surlaie  wavi>s,  ami  i)--e  1 imi  narv  results 
indicate  tfi.it  it  i;.  .i  verv  effective  preproi  essor  ior  I -k  spectra  an.ily.sis 
(IKCUMIi  or  maximum  likelihood  f-k  -.pectra)  ii  it  is  applied  to  all  elements 
of  an  array.  .\|)p  1 ic.at  ion  oi  matciied  filtering  requires  ih.ai  the  signal  wave- 
form bo  known,  wdiich  in  turn  presupposes  knowledge  .if  the  approximate 
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t’l>iienter,  wliii’li  m.iv  he  .icquiied  by  stu)  rt-per  iixl  ik-tt-Lt  ion.  if  tiie  L’picenter 
is  l■nl)v.^,  tiic  rui  I'rd  inp  s ol  .i  lU'iiiby  liirpf  eviiit  (iui  iiL'  usi’ii  .is  tiit*  o;i)ected 
wivoi  liras  . Al  turn.iL  i\'olv,  i!  the  dlsiersion  ibarac  ler  is  t ii’s  ol  tiie  paiii  are 
Kiiiiwii  su:  I icientlv  well,  the  slpii.il  waveloniis  c.iii  he  svnthesized  and  the 
resiiltiiip  wa.  el  on-  used  .is  .1  nsiteheJ  filter. 

oil  .1 1 til  ii.it  i ve  app  1 i I .It  ion  ol  matched  1 ilterinp,  c.'in  he  rei.ative  location 
01  events.  Ii  recordiiiKs  ol  .1  relereiice  event  (prel  eralj]  y ol  an  e>:plosion; 

.ire  .i\ .1 1 1 .ii)I e at  a set  ol  seismic  st.'itions,  the  times  01  ma.'ciin.a  resulting 
i roni  the  matched  I ilterinp  ol  sei.smii  traces  ol  ne.iriiv  events  with  waveforms 
ol  the  rel  ei  em  e event,  can  he  considerei!  a.s  relative  arrival  times  for  the 
purpose  of  event  relocation  in  tin'  penoral  repion  surrounding;  the  reference 
event.  The  technique  also  lias  .1  potent  i.a!  ns  .a  discriminant  since  aElmutlial 
varl.itions  in  the  initial  phases  ol  e.arthquake'S  will  i auso  inconsistencies 
in  the  t i.mes  ol  occurrences  ol  matched  lilter  output  maxima  when  compared  to 
explosions. 

•i. Itched  ! ilterinp  is  essentially  convolution,  and  the  computational 
advantages  of  convolution  or  recursive  filters  on  the  ILhlAC  stated  above 
apply  in  this  c.ise. 

Possible  applications  of  tlie  ILl.lAfi  (Figure  5)  includes  matched  filtering 
01  manv  sites  si.mul  taneousl  v (each  with  a different  matched  filter),  filtering 
several  sets  01  arrav  elements  simultaneously  with  matched  filters  corre- 
sponding to  each  arr.iy,  or  liltering  independent  sites  ysuch  as  LKS.M  sites) 
with  their  own  resjiective  matched  1 liters.  One  can  also  use  matched  filters 
corresponding  to  several  areas  ol  interest  routinely  on  the  data. 

Pill  LTRL 

Ihis  |)rocess  is  designed  lor  a single  tliree-comiionent  set  of  long-period 
data.  It  uses  a nonlinear  weig,hting  .scheme  of  Fourier  siiectral  components 
in  uverlap]iing  time  window.s  to  enhance  Ixjve  or  Rayleigh  particle  motion 
associated  with  a given  arrival  direction  (Simons  1%8).  First  the  three 
components  of  long,  period  recordings  are  rotated  to  obtain  radial  transversal 
and  vertical  motion.  The  rot.ated  traces  .are  broken  up  into  overlapping  time 
windows  and  Fourier  t ran.s  f ormed  , yieldin>;  the  Fourier  coefficients 
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c(i)*cos  27inrT*cl't; 


2 /r 

a (nf)  =-/ 

^ ‘ 0 


9 'P 

b Jnf)  = - / c(T)*si.n  2nnfT(JT; 


where  c(t)  i.s  die  radial,  transverse  or  vertical  component  to  be  analyzed, 

T is  time,  n = 0,1 ,2,3, . . . ,N-1,  Nf  = folding  frequency,  and  f = " fundamental 

liarmonic  of  Fourier  series. 


Using  the  absolute  value  of  a Fourier  component 


A ( n f ) = 2 , 7T  ,2.  ^ 

c a (ni ) + b (nt ) 

c c 


one  computes  three  quantities  used  in  die  weigliting  scheme 


a.)  llie  apparent  liorizontal  azimuth  (the  angle  from  the  radial 


direction) 


3(nf)  = arctan 


A^(nf) 


A^(nf) 


b.)  A measure  of  the  eccentricity  of  the  particle  motion  ellipse 

A ^(nf)  + A^nf) 


'i'(nf)  = arctan 


A (nf) 
z 


c. ) 'ilie  phase  difference  between  the  vertical  and  radial  components 
a(nf)  = G (nf)  - 0 (nf). 


Hie  Fourier  amplitude  coefficients  of  eacli  direction  components  are 
tlien  weiglited  in  the  following  manner 


A'(nf)  = A (nf ) • cos‘*[  3 (nf ) ] • cos  ' 1 'F(nf )- . 21ti  ] • sin  [u(nf)] 
z z 


A'(nf)  =A  (nf)*cos^'[3(nf)]*cos'^t'l'(nf)  - . 2l7r]  • sin*^  [a(nf ) ] 
r r 


A|^(nf)  = A^(nf ) • sin^'l  3(nf  ) ] • s in  ['i'(nf)J’l 


I' 

where  sin  [a(nf)]  0 if  ii  a(nf)  <.  2ii. 
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Tho  A'(nf)  are  the  "weiHlited  amplitude  coefficients".  No  weights  or 
ad  j UK  Lineiu  s are  applied  to  the  pliase  anp.los.  Tlie  exponents  M,  K,  and  N are 
parameters  that  are  read  into  the  program.  Values  of  M,  K,  and  N which  have 
worke<l  reasonably  well  in  practice  range  from  4 to  8.  Note  that  on  the 
vertical  radial  components  all  weigliting  factors  vary  from  1 to  0 as  powers 
ot  sines  and  cosines  depending  upon  tlie  degree  to  which  the  particle  motion 
roseiubles  pure  Love  or  Rayleigh  waves. 

The  cifects  of  the  first  weigliting  factors  (functions  of  f)  are  to 
attenu.ite  transverse  energy  on  tlie  vertical  and  radial  components  and  radial 
energy  to  the  transversed  component. 

Ihe  second  set  of  weigliting  factors  depends  upon  the  angle  f - a measure 
of  tlie  eccentricity  of  the  Kayleigli  orbit  providing  transversal  trace  does 
not  contain  too  much  non-Rayleigh  ty|ie  motion. 

Un  the  vertical  and  radial  traces,  the  angle  desired  (0.21ti)  is  the  one 
corresponding  to  a representative  value  of  the  horixontal/vertical  displace- 
ment ratio  (0.8)  for  fundamental  long-period  Rayleigh  waves. 

The  resulting  Fourier  coefficients  are  subsequently  transformed  back, 
into  the  time  domain  to  yield  transverse  traces  containing  only  Love  motion 
and  radial  and  vertical  traces  with  only  Rayleigh  motion  and  greatly  reduced 
noise  since  the  weighting  scheme  de-emphasixes  noise  which,  even  if  coherent, 
is  liable  to  come  from  a direction  different  from  that  of  the  epicenter. 

The  data  dependent  nature  of  this  algorithm  does  not  lend  itself  well 
to  utilize  the  parallel  computing  feature  of  ILLIAC.  However  if  large  sets 
of  data  need  to  be  processed  each  PF.  can  process  three  components  of  data 
from  a given  location  (Figure  6).  This  may  make  PHILTRE  a practical  pre- 
processor lor  arrays.  Recent  work  by  von  Seggern  and  Sobel  (1974)  indicates 
that  It  is  effective  in  revealing  Rayleigh  waves  hidden  in  noise.  Although 
tiirther  tests  are  needed  to  establish  its  effectiveness  as  a preprocessor  for 
an  f-k  type  detector,  it  utilizes  a neglected  aspect  of  surface  wave  detec- 
t ion. 
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'Jijr'  I ilioml  f-k  Spcctrn 

A m.ixiniuin  likciihiiod  I-k  .spectrum  is  the  mapping  ot  the  power  passed  by 
a .set  01  maxir.mm  Jikelihood  ii  iters  in  the  plane,  A maximum  likelihood  filter 
is  ,'in  optimum  filter  which  is  constrained  to  pass  a plane  wave  in  the  direc- 
tion to  be  li)oked  at  while  rejecting  all  the  re.st  ol  the  ener;;y  present,  in 
a least  mean  square  sense.  It  has  the  mathematical  form  for  a given  frequency 

P(K)  = _i 

u;-l  -t 

li 

.'.here  . is  the  jiower  .spectral  matrix  of  the  sensors,  K is  the  wavenumber  and 
u is  a vector 

i.Kr  iKr  ikr 

u = (e  , e . e “) . 

Ihe  position  vector  of  the  i'th  sensor  of  the  array  is  r., 

Ihe  maximum  likelihood  f-k  spectrum  i.s  one  of  a wide  family  of  high- 
resolution  spectral  estimators.  It  is  characterized  by  reduced  side  lobes 
and  hig.her  resolution  as  manifested  in  the  reduction  of  the  width  of  the  main 
lobe  when  compared  to  the  simple  frequency  domain  beam  used  in  FKCOMB,  The 
processor  requires  the  estimation  of  the  inverse  of  the  input  spectral  matrix; 
there  are  fast  practical  ways  to  make  this  estimate,  after  which  the  multipli- 
cations with  the  various  u vectors  can  be  done  rapidly  by  using  all  64  parallel 
|i roc essor s . Ihe  parallel  teature  can  al.so  be  used  to  Fourier  transform  the 
individual  seismic  trace  segments  simultaneously.  Algorithms  are  available 
to  estimate  the  inverse  ot  the  array  spectral  matrix  without  actually  inverting 
a I'hitrix  (.1.  W.  Woods,  personal  communication,  1972), 

if  the  detection  ol  surface  waves  from  a known  epicenter  is  desired,  the 
range  ol  search  in  the  k plane  is  reduced.  Moreover,  the  absolute  value  of  k 
is  I i.xed  lor  a given  frequency,  since  the  surface  wave  phase  velocity  for  a 
given  I'reciuency  at  a given  array  site  can  be  determined.  Matched  filtering 
or  FHll.IRI-.  can  be  u.sed  as  preprocessors  to  this  processing  scheme  to  utilize 
the  dispersion  and/or  the  particle  motion  characteristics  of  the  signal  and 
reduce  the  false  alarm  rate.  The  most  practical  way  to  use  the  ILLIAC  com- 
puter is  to  ajiply  sixty-four  u vectors  simultaneously  using  the  same  estimate 
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oi  tiu*  coinp.uaMon  hv  a la  lo.  oi  h\  relativi^  to  sequential  processing  ami 
is  the  most  etiicient  lot  the  computation  ot  tinelv  spaced  values  in  the  f-k 
plane  needed  hv  this  h 1 ph- re-u,  | ut  i on  process.  A How  diapram  in  1-ipure  7 shows 
how  the  unique  i)arallel  coniputinp  feature  ol  11. MAC  can  he  used  to  Increase 
the  etiiciencv  of  comput inp  maximum  likelihood  f-k  spectra. 

KKCOMli 

IKCUMB  is  a last  f-k  analysis  i)ioprani  that  was  tirst  used  in  an  automatic 
processinp  svstem  for  mi c roharae raph  data  (Smart  and  I’linn  1971).  it  has 
since  been  adopted  for  use  with  1.1'  seismic  data.  It  computes  and  finds  the 
maximum  of  the  fun.tion 

1 1 

~ ( ■ ) J • exp ( ik*  r ) " 

n=l  " " n 

which  is  ossentiallv  the  power  in  tiie  Irequency  domain  beam.  Here  is  the 
angular  irequency,  A(.)  exp  ( i ( „>  ] is  tlie  Fourier  transiorm  of  the  n'th 
seismic  trace,  N is  the  number  of  traces,  and  exp(ikr^)  are  the  components 
oi  the  vector  u in  the  previous  section.  The  maximum  of  the  function 
can  bo  associated  witli  tlie  presence  oi  a signal.  The  F test  is  used  to 
determine  whetlier  a signal  is  present. 

Ihe  methoils  take  advantage  of  tlie  fact  that  the  signal-tn-noise  ratio 
varies  with  frequency,  so  heamiorming  is  done  frequency  by  frequency.  Also, 
by  staving  in  the  frequency  domain  a groat  many  beams  can  be  examined 
rapid]'.’,  tlie  number  being  limited  on!’,  hv  the  resolution  cell  of  the  array 
response.  The  Jow  resolution  of  tlie  process  is  actually  an  advantage  when 
one  desires  to  search  f-k  space  rapidly,  since  the  wide  main  lobe  of  the 
process  en.ahles  one  to  use  a wide  grid  spacing  in  the  search. 

The  nximutti  and  velocity  of  a signal  need  not  be  assumed:  one  merely 

accepts  the  beam  with  ma.ximun  jiowor.  Ibis  fact  is  important  for  signals  such 
as  long-period  seismic  sun  ace  waves,  which  not  only  are  dispersive  (l.e., 
their  phase  velocities  varv  with  Irequency)  but  whose  arrival  azimuth  may 
also  vary  with  frequency  because  of  lateral  inhomogeneities  in  their  paths. 
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DATE  F.TOM 


SITE  1 SITE  2 SITE  3 SITE  < SITE  5 


PREPTOCESSO.T  (OPTIOMAL) 
(ilVrCHED,  FREQUENCY  FILTER) 
PHILTRE,  AXIS  ROTATIO.I) 


FAST  FOURIER  TRAIISFORM 


COMPUTE  INVERSE 
OF  SPECTRAL  I«\TRIX 
{SINGLE  PE  CPERATIO.I) 


COMPUTE  f-k  SPECTRA  USIIJG 
VARIOUS  BEE  VECTORS  FOR  A 

SET  OF  '.lAVERUriBERS  kj.  k^ 

k^  IN  SETS  OF  64. 


I inure  7.  rowputation  of  Maxinum  Likelihood  f-k 
Spectra, 
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Since  ihc  main  .uivaiU.ij'e  oi  the  l-KCOMh  metluHl  is  ihe  imssibiliiy  of  fast 
search  in  the  vavonumher  vector  s;!ace  at  a ^;ivep  ire(|uenry,  chanjtin^;  1 re- 
cjiiencies  as  the  seatih  retjiiiies,  we  pro)' i aniineiJ  tlie  processor  to  operate  on 
tour  successive  time  windows.  Ihis  uses  the  IIJ.IAC  most  effectively 
for  si(;nal  detection.  ihe  otiier  ivpe  t)i  apj)  1 i cat  i on , searchin);  sixty-tour 
frequency  levels  siriui laneousi v on  the  same  lime  window,  is  not  so  efficient, 
since  not  all  ol  tlie  tiecpiency  ti.iiiJs  may  he  net'de-d  lor  the  search  in  a )^iven 
iterat ion. 

Details  about  tiie  pro)',rainmine  of  IKt.otlii  will  he  piven  in  the  latter 
part  of  this  ro])ort. 
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KKCiMB  AI.(.()KITHM  UKSICN 


' >1' 

11,0  KK(.>MB  can  l,o  di-  Idod  into  iho  following  steps: 

1.  Input  raw  Ion»;  ,,er  lod  data.  Separate  it  by  array.  Extract  the 
Ion,;  period  data  samples  and  the  timin,;  words  associated  with  those  samples. 

2.  l.ivide  the  input  into  time  wind,s.s.  As  oriBlnaUy  input  the  data 
is  ordered  in  tlie  lol  lowing  manner: 

H1.1>.1(2.1i....I(N,1),  l(l,d),i(2.2)....T(N,2l,  1 (i..S>  ,T  (2.  S) , . . .T  (N  ,Si 

wt.ore  I(i,j)  represents  the  data  sample  irom  channel  i at  time  j,  N is  the 
number  of  chamw^ls  and  S is  the  number  of  time  periods.  After  division  into 
time  windows  the  data  is  ordered  as  lollows: 

I (N,rw),l  il,TW*l),l  (1,1W*-2),...T(1,2TW),... 

where  TW  is  the  tine  window  lenKth  and  T(i.j)  and  N are  as  above. 

i.  Ihe  data  is  converted  f rcm,  the  raw  data  format  to  the  internal 
representation  ol  the  machine  used.  (.1  itches  or  spikes  are  removed  and  dead 
<r  noisy  channels  are  detected  and  removed.  (Portions  of  this  step  may  be 
jierJormed  berore  step  2.) 

*.  A Tourier  transform  is  applied  to  each  time  window.  Alter  FR  the 
‘..it.i  is  -if a.s  lollows; 

Ml.1,1  l.Kl.l  .-’>.••.‘•■(1  ,l,rWl,K(l,2,li.F(l  .2,2),.  ..F(1,2,TW),... 
m.;;.l).F(l.N.2)....F(l,N.lW,.F(2.1.1).F(2.1,2)....F(2.1.TW). 

K2 , 2 . 1 ) , F(2,2. 2) . . . . F (2,2  ,iW)  .... 

w.ure  FU.i.k),  the  Fourier  transform  output,  represents  frequency  k,  channel 
. t i.'.e  window  i. 

i.  Ke-order  tta.-  data  so  that  it  is  arranged  by  trequency.  it  is  then 
arranged  as  lollows: 

K1,1,1),K1,2,1),...F(1,N,1),F(1,1,2),F(1,2,2),...F(1,N,2),... 

F(1  ,1 ,TW),F(1,2,TW) ,...F(l,N.rW) ,F (2, 1 , 1) , . . . 


where  k(i,i,k),  IW,  and  N are  defined  as  above. 


b.  Search  frequency  - wavenumber  space  lor  power  maxima. 

Data  KdltlnK  Module  One  (IIKMI) 

Since  step  one  is  a process  common  to  all  seismoloKical  alKorithms 
and  because  lar«e  input/output  buffers  are  required  it  was  coded  as  a stand 
alone  module.  The  input  to  this  module  (DKMl)  is  the  raw  data  as  read  from 
a low  rate  tape.  The  format  of  input  records  is  shown  in  Figure  8.  The 
output  consists  of  several  files,  one  per  array,  containing  only  the  data 
samples  applicable  to  lon^  period  processing.  The  data  movement  required  to 
isolate  and  properly  structure  this  data  is  nonparallel.  There  are  no  general 
structures  repeated  often  enough  to  allow  efficient  use  of  the  ROUTE  instruction, 
The  CU  is  used  to  move  one  word  at  a time  between  input  buffers  and  output 
buffers.  (Actually  the  HIN  instruction  is  used  to  move  blocks  of  eight  words 
to  the  CU.)  A description  of  each  array  format  is  given  in  the  block  data 
subroutine  initialization  ol'  the  vector  (NTRL. 

The  reorderlne  of  date  In  „e,.s  2 and  5 are  not  reqnired  If  all  data  la 
available  on  a randot,  acceaa  device,  Unce  It  reflects  the  order  In  .hlch  the 
data  „UI  be  accessed.  It  Is  necessatv  on  1I.UAC  since  the  sl.te  ol  core  and 
long  disk  access  time  prohibit  random  access. 

Assuming  approximately  20  channels  for  each  of  three  arrays,  each 
sampling  at  the  rate  of  once  per  second,  one  twenty  four  hour  tape  contains: 

3 arrays/sample  * 20  cliannels/array  * 1 sample/second  * 

60  seconds/minute  * 60  minutes/hour  * 2A  hours/day  = 5,184,000  channels. 

Moving  each  sample  involves  two  memory  accesses  (one  load  and  one  store). 

Given  that  a memory  access  from  the  CU  takes  approximately  .5  microseconds, 
the  total  time  spent  in  memory  accesses  by  DEMI  to  process  twenty  four  hours 
of  data  is  on  the  order  of  5 seconds.  This  is  small  enough  that  more  com- 
plicated algorithms  which  may  have  permitted  use  of  the  ROUTE  instruction 
were  not  considered. 

Twenty  four  hours  of  d.ita  is  approxlnately  ten  to  the  eight  bits.  In 
order  to  read  these  Into  core  without  losing  a great  deal  of  time  waiting 
for  disk  access  a buffer  of  128  rows  (512,000  bits)  of  core  Is  used.  200 
disk  accesses  are  required  for  Input.  This  takes  up  to  eight  seconds.  Since 
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there  are  several  oiitinit  tiles,  the  outi)iii  bul'fers  are  somewhat  smaller.  The 
total  size  of  the  output  is  smaller,  since  at  least  half  of  the  input  is  not 
used  in  long  period  processing.  The  i/O  time  spent  in  oiitimt  is  therefore 
approximately  equal  to  that  spent  in  input  even  though  the  output  buffers 
are  smaller. 

The  actual  movement  of  data  hv  DMSl  is  done  within  three  nested  loops. 

The  outermost  loop  is  gone  thru  once  for  each  input  record. 

The  next  inner  loop  is  gone  through  once  for  each  time  scan  in  each 
record.  The  innermost  loop  is  gone  thru  once  for  each  channel  per  time  scan. 
All  buffering  is  handled  by  an  input  routine  and  an  output  routine  called 
once  for  each  channel  to  transfer.  In  order  to  reduce  overhead  spent  in 
subroutine  calls  it  may  be  necessary  to  recode  calls  on  these  routines  as  in 
line  code. 

DfcMl  transposes  data  in  a serial  fashion.  It  is  coded  so  as  to  minimize 
time  lost  in  disk  and  memory  accesses.  It  puts  array  data  in  a standard 
format  to  reduce  the  size  of  the  data  and  allow  the  straight  forward  execu- 
tion of  subsequent  modules. 

Data  Editing  Module  Two  (DEM2) 

Steps  2 through  5 are  performed  by  i)IM2.  The  primary  reason  that  this 
module  was  coded  separately  from  step  6 was  to  shorten  coding  and  debugging 
time.  The  relatively  small  amount  of  core  memory  available  in  each  PE  would 
necessitate  the  overlaying  of  various  vectors  used  by  step  6 and  those 
included  in  DEM2  if  all  were  included  in  one  module.  The  1/0  times  spent 
writing  the  output  from  DLM2  and  reading  it  in  before  step  6 would  be  saved, 
but  this  time  is  estimated  to  be  less  than  5 seconds. 

Steps  2 through  5 are  performed  one  time  window  at  a time.  A complete 
multi-channel  time  window  is  taken  through  steps  2 through  5 and  the  resulting 
output  placed  in  an  output  buffer  before  the  next  time  window  is  processed. 

One  multi-channel  time  window  consists  of  approximately  20  channels  of 
up  to  512  samples  each  or  approximately  10,000  data  items.  During  step  2 it 
is  impossible  to  Include  a complete  multi-channel  time  window  within  one 
processing  element  memory.  It  is  possible  to  include  a single  channel  time 
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window  wiihin  one  piitc  css  i iik  clement  mejiiorv,  luit  cine  to  the  variable  number 
ol  iluinnels  used  lot  e.uh  time  window,  kee|)in);  track  oi  which  channels  and 
time  windows  have  been  processed  i.s  complicated,  thoii)’h  feasible. 

An  alternate  .ippmach  is  to  use  64  l‘i 's  to  process  each  time  window. 

An  1-l-T  routine  is  available  (written  at  tiu'  I'niversity  of  Illinois;  which 
utilizes  all  ()4  IT's  to  perform  one  H-T  which  runs  very  close  to  64  times 
taster  than  one  l‘i;  would  ilo.  Conversion  to  float  inp,  point  involves  no 
interaction  between  processinn  elements.  Dep 1 i tchinp  involves  the  comparison 
of  eacli  sample  wath  the  previous  and  ne:-:t  s.Enple.  With  this  data  arrangement 
these  samples  are  in  adjacent  PK  ’ s and  the  HOl'Ti;  instruction  can  be  effectively 
utilized.  I he  oripinal  structurinp  of  the  data  into  time  windows  (step  2) 
and  the  final  transposition  (steji  3)  are  each  |>erformed  serially  by  tlie  CL!, 
so  are  not  affected  hv  the  data  arraiu;emeni  chosen.  Spread  inp  time  ...'ndows 
across  I'L'.s  was  the  approach  chosen  for  steps  2 throu);h  3. 

Step  2 thus  consists  of  estractinp  ilminp  information  from  the  input  and, 
usinK  this  information,  form  t ime  wi  lulows.  I.ach  time  window  is  spread  across 
the  I’K  s,  occupyinp  between  one  and  cipiit  rows  per  channel,  de|)endins  upon 
the  time  window  size  in  use,  Overlappinp  of  time  windows  is  performed  by 
retaininp  whatever  part  of  the  most  recent  time  window  is  still  of  interest 
and  usinp,  the  ROl'i  K instruction  to  hack  it  up  properly.  Tor  this  reason, 
the  buffer  in  whicli  time  windows  are  built  Ls  alternated  between  two  halves 
of  an  array  so  tluct  tiic  last  time  windcjw*  built  is  not  overwritten. 

Conversion  to  internal  floatinj;  point  format  is  the  first  step  per- 
formed once  tlic  time  windows  have  been  formed.  bach  iT  converts  all  samples 
within  its  memory  and  no  inter  I'K  communication  is  required.  Deglitching 
is  performed  hv  routing  the  values  from  adjacent  IT.'s  and  adjusting  them  if 
a glitcli  is  encountered.  (See  ])rogram  documentation  for  exact  procedure.) 

The  rowsum  procedure  detcrihed  in  section  2 is  used  in  the  variance  calcula- 
tion, since  a summation  across  IT.’s  is  required.  The  FFT  is  then  performed 
and  the  frequencies  prepared  for  output.  Duo  to  the  fact  that  not  all  fre- 
quencies output  by  FFT  are  of  sei  smo  log  leal  interest,  the  output  from  step  4 
is  much  smaller  than  the  input  to  step  4.  This  data  reduction  is  significant 
in  that  after  FFT  a multi-channel  time  window  consists  of  approximately  20 
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channels  and  20-JO  I reqnencies  and  will  111  wiihin  one  prcjcessinK  element 
memory.  I’lacinji  each  separate  time  window  whollv  within  one  I’r.M  is  very 
convenient  since  the  search  of  Irequencv  wavemunher  space  lor  one  time 
window  is  completely  sell  contained  .and  independent  of  other  processin;;  (see 
Figure  9).  In  step  5,  the  l-FT  output  is  written  into  one  I’l:'!  ol  an  output 
core  buffer.  This  is  done  serially  by  the  CL'.  When  the  butler  is  full  it 
is  written  to  disk. 

FKCOMB  Algorithm 

Since  each  Lh  completely  cont.ains  one  lime  window  af  ter  slej)  5,  the 
algorithm  used  in  the  search  of  frequency  wavenumber  space  is  essentially  the 
same  as  that  used  in  the  serial  version.  The  program  reads  the  data  file 
created  by  DKM2,  which  has  been  arranged  as  shown  in  Figure  9.  Kach  time 
window  contains  the  frequencies  ot  interest  and  the  algorithm  is  executed 
in  parallel  on  the  data.  A search  for  maximum  power  is  made  on  a coarse 
grid  .■’.nd  then  a series  of  finer  grids  is  se-arched  simultaneously  in  all 
processing  elements  until  a maximum  is  found.  In  a given  PK  the  mode  for 
that  PE  is  disabled  until  a maximum  is  found  in  all  other  PH's.  The 
Fisher  statistic,  period,  s'lgnal  azimuth  and  velocity,  and  associated  para- 
meters are  calculated  and  stored,  and  the  process  is  continued  on  the  next 
time  window  of  data.  The  design  of  the  algorithm  w.'is  straightforward,  and 
the  reader  is  referred  to  the  prograi:;  documentation  (Kerr  and  Wagenbreth, 
1974)  for  a more  detailed  discussion  of  the  software. 
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The  operational  aspects  of  usin^  ]U,JAC  differ  si;;nif icantly  from  those 
of  otlier  machines.  in  addition  to  die  parallel  arcliitecture,  Lliere  are  two 
other  characteristics  whicli  are  important  considerations  for  the  user  of  tlie 
ILLIAC  system.  First,  all  of  the  support  software  sucli  ns  editors  and  com- 
pilers run  on  processors  other  than  I1,LIAC.  lliere  is  currently  support  soft- 
ware available  on  DEC,  IBM,  and  Burrouglis  machines.  The  clioice  of  which 
machine  and  software  to  use  is  an  integral  part  of  system  development,  for 
ILLIAC  is  accessible  only  via  the  ARl’A  Network  and  is  routinely  used  remotely. 
Tlie  bandwidth,  availability,  and  reliability  of  tlie  network  directly  affect 
the  performance  of  the  ILLIAC  system  as  seen  by  a user. 

Program  Entry  and  Storage 

A basic  requirement  for  any  long-term  coding  effort  is  a reliable  file 
system  permitting  easy  access  and  modification  of  source  codes.  Two  basic 
options  were  available  in  using  the  ILLIAC  system.  One,  used  by  several 
ILLIAC  coding  efforts,  is  to  maintain  files  on  a host  computer  and  transfer 
the  files  to  the  ILLIAC  system  via  the  ARPA  Network  whenever  necessary.  Tlie 
second  is  to  utilize  the  Tenex  file  system  and  editors  included  in  the  ILLIAC 
system.  Tlie  first  api  reach  required  frequent  ARPA  Network  transfers  and  a 
reliable  and  economical  host  machine.  Since  such  a host  was  not  available 
to  SDAC,  the  Tenex  file  system  was  used,  and  was  found  reliable  and  convenient. 
No  work  was  lost  due  to  disk  or  file  system  failures  during  tlie  duration  of 
this  project.  Tlie  editor  DED  fulfilled  all  requirements  regarding  both 
modification  and  examination  of  source  files. 

Languages 

Three  languages  are  available  for  preparation  of  ILLIAC  code.  Tliere  are 
two  high-level  languages,  GLYPNIR  and  CFl),  and  an  assembly  language,  ASK. 

The  large  amount  of  coding  necessary  made  the  use  of  assembly  language 
impractical  except  for  specific  portions  where  bit  manipulation  or  efficiency 
made  it  a necessity.  Tlie  majority  nf  coding  was  done  in  high-level  language. 
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A comparison  ot  iho  syntax  and  somatulcs  oi  CI.M-.mk  and  ( KU  lovt-aU-d  u,.- 
lollowinn  Hit-ni  I icant  d i f 1 oreiico.s : 

1.  i.aso  OI  nndorstandinn  - CKI)  looks  nuicl.  like  loKlKAN  and  is  easily 
inurpreud  or  learned  by  a sciontiii,  prop.ranmar.  (d.Yl'MK  resembies  AHol. 
..nd  is  somewhat  more  .onlnsinn  and  dilliciilt  to  learn. 

d.  base  OI  codinp  - Once  learned.  Cl.Yi'MH  permits  laster  and  dearer 
■odiin.  than  CKI).  Cl.iPMKVs  macro  lacilities  are  a convenience  not  provided 
bv  (Kb.  (d.Yl'NIl;  has  some  bibber  level  ronstrnct.  which  require  several  (dll 
statements  to  implement. 

t.  Kllidency  - CKI)  produces  more  etii.  ient  code  than  (d.Yi’Ml;  in  manv 
instam  es. 

iiie  two  lannaaites  are  very  similar  in  ti.eir  treatment  ol  unique  11.;,1A( 
characteristics  and  both  provide  all  lacilities  necessarv  lor  il.e  implementa- 
tion ot  seismic  analysis  programs.  Certain  types  oi  code  are  !,et ter  suited 

to  one  ian«ua,;e  chan  the  .tber,  but  consideration  o!  tire  syntax  and  semanti 

alone  indicated  no  clear  prelerence. 

n.e  choice  ol  lan^uap.e  ultimately  depended  upon  the  support  and  aval  I- 
a.  ility  Ol  cdi'Mk  and  (.Ki).  Cl.YKMK  is  supported  bv  the  Institute  tor 
■-  ained  Computation  as  part  ol  the  II  I 1A(  svster.  It  is  implemented  or  a 
■wn.ucn:,  c,  .0  Kuated  at  the  lid.iAC  computer  .enter  and  must  be  ao  esse.t  ... 

the  lU.  V ,;tc!,  que.ie  ias  discis.ed  b.  low,.  (b.,  i,  implemc-nted  on  the 

lo.oted  at  riASA  Ames  ke  .e.ircii  . eni.r.  it  i...  .ocessiblc  roittn.  l 
la  the  AKK.-.  hetwork.  n.e  source  t..r  ( KI,  is  t r.insportahle  and  a version  ..1 
■ is  available  on  the  ftl.A  II;m  IM./d.  Also  s.q.porte.l  o„  the  .Mtes  IBM  Ibn,.,. 
i a tr.mslat..r.  t Kl,>. . which  translates  crq,  ,o  Kortr.m.  With  s..ne  mod  i ; i c i- 
t.on  due  to  !/„  ditferonces  ..nd  inscrt...l  asser-.I.l-.  Ian,uia,te  code,  (.KI,  provr-ims 
mav  be  tr.in.l.ated  to  I Hit  lortr.in.  Ibe  translator  is  desired  to  yeneiMte 
-cle  equivalent  t..  that  generated  bv  (.i  |,  ,„r  ii|.lA(.  iToyrarrs  ,r.av  tt.en  b. 
Jebu.q.ed  andlste.l  on  an  Ih.M  IM,  rather  than  ..n  ll.l.l.(t.  ,,.„x  is  not  d.sicned 
to  replace  ll.l.l.A(  in  producti..n  ru.de  sine  tl.c  Knki  KA.N  penerated  by  writing 
a CKI,  pr.n,.r.im  ,.nd  translati.n  it  will  n..t  he  nea,  Jv  as  et  t i.  ient  as  .odim; 
in  b',l<H;A;;  .lire,  t Iv. 
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Due  to  the  superior  avaUahility  of  LTD  and  the  existence  of  the  CFDX 
translator  the  decision  was  made  to  implement  the  I'K(  OMli  algoritlim  in  CFD. 

As  actually  experienced,  tlie  availability  of  CFD  was  not  as  good  as  had 
been  hoped,  for  several  reasons;  First,  availability  of  the  Ames  115M  360  is 
very  poor.  Between  the  hours  of  8:00  am  and  12:00  midnight  (I’ST)  use  of  the 
machine  by  non-priority  accounts  is  restricted.  During  the  eiglit  remaining 
hours,  the  requirement  that  botli  tlie  ii.UAC  Tenex  system  and  the  Ames  IBM 

J60  be  operational  for  file  transfers  caused  mucii  lost  programmer  and  com- 

puter time.  The  hours  were  also  inconvenient.  Tlie  L'CIA  version  of  CFD,  due 
to  lack  of  overlays,  requires  400K  core  and  runs  in  a slow  queue  (6-8  hours 
turnaround).  Efforts  to  Implement  CFD  on  the  8DAC  IBM  J60/44  were  frustrated 
due  to  incompatability  of  the  operating  systems  of  tiie  1 15M  360/44  and  the 

IBM  360/67.  The  large  core  requirement  also  posed  a serious  problem.  it 

was  found  that  the  effort  required  to  Implement  CFD  at  SDAC  would  not  be 
worth  the  convenience  of  an  in-house  compiler.  The  availability  of  ILLIAC 
was  sufficient  (see  below)  to  make  tlie  use  of  the  CFDX  translator  uneconomical 
due  to  the  alterations  necessary  to  accommodate  1/0  differences  and  inserted 
assembly  language  statements. 

Run  Procedures 

Coding  of  the  FKCOMB  algorithm  in  CFD  began  in  April  1974.  ITliat  follows 
is  the  set  of  procedures  developed  for  the  day  to  day  process  of  running  and 
debugging  an  IILIAC  program,  along  with  experience  gained  and  observations 
made  during  the  use  of  these  procedures. 

The  primary  site  at  wtiich  compiles  were  done  was  the  i\mes  IBM  360/67, 

A CFD  restriction  is  tliat  all  subroutines  must  be  separately  compiled.  Our 
code  was  divided  into  three  programs,  each  consisting  of  a main  driver  and 
four  to  six  subroutines.  Initially  all  subroutines  had  to  be  compiled,  but 
thereafter  only  tiiose  with  code  modifications  required  compilation.  Com- 
piling a module  consists  of  four  steps.  First,  after  having  logged  in  on 
the  /\mes  360/67,  the  source  file  is  transferred  over  tlie  ARi’A  network  from 
the  14-Tenex  File  System,  where  the  source  files  are  maintained,  to  the 
Ames  360.  This  process  is  done  interactively  and  typically  takes  one  to  ten 


-35- 


minutt-s  ol  ro.il  time.  depondiiiK  upon  the  lenKlh  ol  the  source  and  the  Joad 
uveruKc  on  end,  ntuliine.  Approximately  one  out  ol  tl.ree  transfers  terminated 
ahnornuilly  and  had  to  ho  reinitiated.  i1,e  failure  rate  increased  Kreatly 
wlien  the  load  on  either  machine  was  heavy.  •II, e next  step  is  to  initiate  tlie 
(.11)  compiler.  ilie  time  between  the  submission  of  a comiiile  and  its  comple- 
tion varied  from  five  minutes  to  several  hours,  ap.ain  dependent  upon  the 
macliine  load.  After  termination  of  tl,e  compile  the  listinj-  generated  by  tlie 
CFl)  compiler  is  examined  witl,  the  TSS  editor,  KliDlT,  to  clieck  for  syntax 
errors  or  other  abnormal  termination.  If  errors  are  detected,  tl.ey  are  cor- 
rected (heinn  caret ul  to  make  the  same  corrections  to  the  original  source  at 
14-Tenex)  and  the  compile  reinitiated.  After  a successful  compile,  the  /VSK 
assembly  language  source  module  is  copied  back  to  I4-’Ienex  via  network 
transler.  This  file  is  usually  several  times  larger  than  tl.e  original  source 
and  the  time  taken  to  trausier  the  file  is  several  times  longer  than  that  for 
the  source.  If  several  subroutines  are  to  be  recompiled,  tliis  process  can 
consume  several  liours.  When  only  small  changes  are  necessary,  tliis  time  can 
be  saved  by  changing  tne  assembly  language  code  directly  witl,  the  text  editor 

nt  U-Tenex,  again  being  careful  to  rmake  the  same  changes  to  the  original 
source. 

Once  tlie  necessary  assembly  language  modules  have  been  created,  a batch 
job  is  submitted  at  14-Tenex  to  perform  the  following  tasks: 

1.  Assemble  the  /VSK  modules 

2.  Linkedit  the  resultant  relocatable  modules 

3.  Create  a disk  map  file  describing  the  actual  layout  of  any  ILLIAC 
disk  areas  to  he  used  hy  tliis  run 

4.  Allocate  the  rmap  file  created  in  the  last  step 

5.  dove  any  input  files  rocpiired  to  the  appropriate  IhhlAC  disk  area 

h.  Klin  the  11.1.1  AC  code 

7.  .Move  any  output  from  the  appropriate  11.1, 1,\C  disk  area  to  the 
14-Teiiex  file  system 

d.  Release  the  ll.Ll,\C  disk  areas  used. 
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Turnaround  for  il.iJAC  liau;h  jobs  improved  1 icanL  1 y I roni  April  t.o 

June  1974  , but.  was  tilways  sul)ject.  Lo  unpred  iclab  le  riiRt.uat.ions  and  delays. 

In  April  no  more  tluin  Lliree  or  four  turnarounds  per  week  could  be  expected. 

Ky  June  tills  had  increased  to  two  per  day  if  submitted  takinj;  into  account 
the  hours  scheduled  for  hatch  jobs  to  be  run.  This  required  .ilmost  constant 
monitoring  of  the  batch  queue  between  tlie  liours  of  9 am  and  9 iim/KDT. 

Turnaround  was  significantly  imiiroved  during  the  montli  of  June  by 
relocating  to  California  and  working  at  the  Ibl.lAC  .site.  Considerable  effort 
was  made  by  the  Ibl.lAC  user  support  group  during  tliis  time  to  insure  that 
SDAC  jobs  were  g,iven  priority,  and  considerable  progrtss  was  made  during  this 
per iod . 

Analysis  of  the  results  of  an  111 lAC  run  is  possible  by  three  methods. 

The  primary  and  bv  far  most  convenient  means  is  an  unsophis  ticateil  form  of 
formatted  I/O  called  "DISIM.AV".  The  output  is  readable,  and  predicted 
answers  can  be  checked  against  actual  results  with  its  use.  In  the  case  of 
unexpected  results  wliere  suitable  displays  were  not  generated  to  provide  clues 
to  the  source  of  the  error,  we  found  it  necessary  to  examine  the  dumj)  files. 
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CONCl.US  1 ONS  AND  UKCOMMl'NDAT  1 ONS 


Sc-lsinlc  1‘ ruccss  i 11^;  on  1 U.  1 AC, 

'Hit-r  Il.LlAi;  cnnipiUfr  proj; rammed  to  jierform  seismic  processiiij;  on  lar);e 
data  i)a.ses  can  be  a v.ilualilt*  tool  in  the  ii<>vei oi)menl  of  seismic  event  detec- 
tion and  discrimination  iH'ocedures.  It  is  feasil)le  to  imi^lement  some  existint’ 
alnoritiims  on  tiie  ll.l.lACi  wiiicli  are  not  currently  used  to  process  lar^e  data 
liases,  or  some  algoritiims  wiiicli  are  jiroposed  hut  not  tested  due  to  a lack  of 
computing  power.  Our  exiieri  ence  with  one  algorithm  (I'KCOMB)  wliich  is  repre- 
sentative of  seismic  analysis  prof;rams  sliows  tiiat  a major  benefit  of  tiie 
Ihl.IAC  to  seismic  processinj;  Is  its  aliility  to  operate  in  parallel  on  sixty- 
four  different  data  streams,  thereliy  reducinp,  tlie  time  required  to  process 
lar^e  data  bases.  Efficiently  arranj^inp,  these  data  streams  for  tiie  processing; 
element  memories  is  an  important  consideration  for  desijtninR  any  seismic 
algoritiims  for  tlie  11.1,1  AC. 

It  is  feasible  to  iiroyr.am  II.I.IAC  to  perform  the  algorithms  reviewed  in 
tliis  study:  convolution-recursive  filtering,  PHILTRE,  matched  filtering, 

heamforming,  and  maximum  likelihood  f-k  estimation.  Since  a major  factor  in 
jirogramming  any  of  these  ;dgor  i thins  is  the  data  arrangment  in  core,  a more 
detailed  study  of  the  data  configurations  for  these  algorithms  would  be 
needed  to  optimise  the  use  of  the  computing  power  of  I LEI  AC.  One  algorithm 
(EKCOMll)  was  studied  in  detail  iuui  implemt-nled  on  ILLIAC  IV.  Data  editing 
schemes  were  devised  for  I'KCUMIJ  which  can  he  usi-d  with  appropriate  modifica- 
tions for  .ill  the  sei  smolog  leal  al);orithm.s  we  reviewi-d. 

'IVo  imle jieiident  use.s  for  ILl.I.AC  are  suggesteii.  First,  FKCUMIJ  and  other 
algorithms  now  used  selectively  could  be  run  routinely  on  larger  data  bases 
to  better  provide  the  services  they  already  g,ive  on  conventional  nuichines. 
Second,  experimental  methods  impractical  to  lest  via  conventional  maciiines 
could  lie  tested  on  11, Li  AC.  'llie  ex])erience  of  implementing  FKCOMB  illustrates 
that  the  desig.n  and  coding  of  new  algorithms  for  ILLIAC  is  not  significantly 
more  difficult  than  for  serial  machines.  'llie  only  phase  not  exjier  imentally 
explored  by  this  effort  are  the  opi'Cational  problems  of  manipulating  the 
large  amounts  of  liata  involved  in  routine  processinj;  of  long  and  short  period 
data  on  ILLIAC. 
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!v)  1'i.ix  inii  ill'  f 1 1 i r i I lu".  , i!k-  ttni-  consiinii'tl  oxt’i'u  L i n>;  an.'ilysis  al)'.niiiir 
sluHilil  hi'  si'.'.nil  ii'iinilv  -if, iter  ih.iii  (In,-  i i im,'  ri‘f|ii  i I'l'il  lor  data  (.“d  i L i ir,.' . 
coni)  i iKi  t i on  oi  a 1 c.o  r i t Inns  such  as  n.atcliod  lilierin)'  and  I’KCOMli  ri>(|uiro  an 
order  oi  nay, nilude  inori  |ii  oci’Ss  in.i’,  timo  on  a y.iven  memory  load  than  l-'KCo.'h; 
alone,  .iml  would  therel)'/  util  life  ll,I.IA(,'  r»ire  e I I i c ien  1 1 y . 

I'  roi'.r:  inniiiy,  1'  Kj.  X )!■  Ijl 

The  I'ollowinit  points  represom  our  1 indin};:'  in  dovelopinj;  I'KCOMTi  sot't'.-.ai. 
on  the  ll.l.lAC. 

1.  Fa-ter  .ind  more  reliable  netvvork.  tile  transfer  between  1 4-Tene.y 

;ind  otlu“r  nosts  such  as  I'Cl.A,  Ames,  I'SS  and  S1)A(!  would  e>;|iedite  the  prop  r.ii'!:  in 
and  use  of  tile  ILl.l  iC  systi".)!. 

2.  'Iliere  is  no  clear  preference  betw»*en  Cl'l)  and  flLVl’NIK  indicated  by 

our  e>:pe  r i ences . Pie  possibility  of  implement  ini;  ■•t-  bUAC  or  upi;radini; 

servile  at  L'Cl.A  should  be  investigated,  and  .in  experiment  made  in  the  use 
o!  (ILYl'MK  iieiore  any  lony  ranp.e  decision  is  made, 

1.  Soitware  debupp.in};  aids  presently  av.iil.ible  for  ILbLAC  pro|;rams  arc 
minimal.  Additional  debnypini;  .aids  would  lessen  the  tusk  f)f  ILLIAC  pronr.aintai  n 
User  implementation  of  such  aids  on  tlie  SDAC  host  is  feasible  thouph  .it  the 
cost  OI  considerable  el  fort. 

,','e  estim.ite  that  the  time  .and  el  fort  required  to  desii;n  and  coik 
an  M.I.IAC  propram  is  ro  imre  than  twice  tliat  retpiired  for  a conventional 
mac.hine.  line  to  the  fact  that  1 1.L  I AC  is  not  fully  operational  at  present 
.and  tin-  necess  i t v to  handle  tiie  l.arpe  amounts  of  data  inherent  in  the  use 
oi  an  till  r.il. 1st  rtachine,  the  time  and  effort  reipiired  currently  to  debup,  .an 
11.1.1  AC  jirovr.im  m.av  he  as  much  .is  four  times  th;U  recpiired  for  a convention.!! 
laarh  i lie . 

b.  Routine  jirocess  i up  oi  2h  hours  of  lonp-period  seismic  dat.i  is  not 
feasible  at  present  due  to  the  restricted  availability  of  the  ILI.IAC  p roces si ■ r 
and  tlu'  inability  at  tliC'  system  to  handle  the  larpc  amounts  of  data  efficient  1 
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Cl.O.SSAKY  ()|-  TKKMS 


Ali'OI.  - lliKlicT  Jovel  alnuritini  I !m^<,ua^■e  on  serial  computers. 

M.I’A  - Alask.in  I-onp  I’eriod  Array. 

ASK  - Assemliiy  IanKua);e  tor  Il.l.IAC. 

a.sSFMRI.I.  - (tonvert  1 rom  mnemonic  code  to  machine  code. 

at  I - lilts  per  inch,  tlie  unit  for  measuring;  inaKnetic  tape  recording  density. 
(ID  - A Fortian-1 ike  higher  level  language  developed  by  the  Computational 
Fluid  Dynamics  group  at  NASA/Ames. 

CFliv  - A progr;im  wliicli  translates  C.Tl)  code  into  Fortran  code. 

Dl.t.  - Digital  Fquipmeiit  Corporation. 

DFI)  - The  text  editing,  subsystem  on  IA-iene>;. 

DI.Ml  - Data  Kditing  Module  1,  the  first  data  editing  program  module. 

DI.N2  - Data  Kditing  Module  2,  the  second  data  editing  program  module. 

DISI'KAY  - Macros  whicli  invoke  subroutines  permitting  print  output  of 
intermediate  results. 

DL'MF  - hexadecimal  representation  of  the  contents  of  core  and  registers  at 
termination  of  program  execution. 

KDI  - Kastern  Daylight  Time. 


MI’  - File  Transfer  I’rotocol.  A protocol  for  file  transfer  from  one  ARPA 
Network  computer  to  another. 

M l - las':  lourior  Transform,  a Fourier  transformation  program. 

FFCOMD  - irequency-k  wavenumber  combination,  the  algorithm  implemented  on 
ll-I.IAf.  discussed  in  this  report. 
l■OHIK/^^  - higher  level  language  on  serial  computers. 

(d.il..  IR  - I lie  higher  level  language  liniil  ewented  for  the  lU.lAC. 

^ Fe r iiii t i o na  1 liusiness  Machines  Corporation. 


Il.l.IAC  IV  - The  parallel  processor. 

IFI.IAC  .Systtsii  - The  complete  computer  system  consisting  of  the  parallel 
proce.ssor,  I’DD-IO's,  Bb70f),  IJNICON,  several  I’DP-ll's,  and  their 
soltware  operating  system. 

M lene>.  - Design. ites  the  I’DP-IO  of  tlie  Il.l.IAC  system,  which  is  the  ARl’A 
Network  host  for  the  Il.l.IAC  system.  14  designates  the  system  with 
respect  to  the  ARI’A  Network. 
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I.ASA  - Large  Aperture  Seismic  Array  in  Montana. 

UNKEDIT  - Merge  several  machine-readable  subroutines  into  one  execulal 
program. 

NASA/Ames  - (also  M,es)  National  Aeronautical  and  Space  Administration. 

Ames  Research  facilitv. 

NET  - Network  Event  Processor. 

NOKSAk  - Norwegian  Seismic  Array. 

PE  - Processing  Element  of  the  lEElAC,  64  of  which  operate  in  parallel . 
PHILTRE  - An  algorithm  using  a non-linear  weighting  scheme  of  Fourier 
spectral  components  to  enliance  Love  or  Rayleigh  waves. 

PST  - Pacific  Standard  Time. 

REDIT  - The  text  editor  for  TSS. 

SDAC  - Seismic  Data  Analysis  Center. 

TENEX  - The  operating  system  for  the  PDP-IO. 

TSS  - Time  sharing  system  for  IBM  360/67. 

UCLA  - University  of  California  at  l.os  Angeles. 

UNICON  - Laser  mass  storage  system  developed  by  Precision  Instruments. 
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